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ABSTRACT 

Inertial waves governed by Coriolis forces may play an important role in sev- 

^ ! 

04 ' eral astrophysical settings, such as eg. tidal interactions, which may occur in 

fT^ . extrasolar planetary systems and close binary systems, or in rotating com- 
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pact objects emitting gravitational waves. Additionally, they are of interest 



C^ ' in other research fields, eg. in geophysics. 



However, their analysis is complicated by the fact that in the inviscid case 
the normal mode spectrum is either everywhere dense or continuous in any 
p% ' frequency interval contained within the inertial range. Moreover, the equations 

Cu ^ governing the corresponding eigenproblem are, in general, non-separable. 

In this paper we develop a consistent WKBJ formalism, together with 
a formal first order perturbation theory for calculating the properties of the 
normal modes of a uniformly rotating coreless body (modelled as a poly- 
trope and referred hereafter to as a planet) under the assumption of a spheri- 
cally symmetric structure. The eigenfrequencies, spatial form of the associated 
eigenfunctions and other properties we obtained analytically using the WKBJ 
eigenfunctions are in good agreement with corresponding results obtained by 
numerical means for a variety of planet models even for global modes with 
a large scale distribution of perturbed quantities. This indicates that even 
though they are embedded in a dense spectrum, such modes can be identified 
and followed as model parameters changed and that first order perturbation 
theory can be applied. 
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This is used to estimate corrections to the eigenfrequencies as a conse- 
quence of the anelastic approximation, which we argue here to be smaU when 
the rotation frequency is sniaU. These are compared with simulation results 
in an accompanying paper with a good agreement between theoretical and 
numerical results. 

The results reported here may provide a basis of theoretical investigations 
of inertial waves in many astrophysical and other applications, where a rotat- 
ing body can be modelled as a uniformly rotating barotropic object, for which 
the density has, close to its surface, an approximately power law dependence 
on distance from the surface. 

Key words: hydrodynamics; stars: oscillations, binaries, rotation; planetary 
systems: formation 

1 INTRODUCTION 

In astrophysical applications inertial waves that can exist in rotating bodies may be excited 
by several different physical mechanisms, most notably through tidal perturbation by a 
companion (eg. Papaloizou & Pringle 1981, hereafter PP) or in the case of compact objects 
through secular instability arising through gravitational wave losses (eg. Chandrasekhar 
1970, Friedman & Schutz 1978, Andersson 1998, Friedman & Morsink 1998). They also 
can play a role in other physical systems. For example, they can also be excited by several 
mechanisms in the Earth's fluid core with possible detection being announced (Aldridge & 
Lumb 1987). 

For rotating planets and stars that have a barotropic equation of state these wave modes 
are governed by Coriolis forces and so have oscillation periods that are comparable to the 
rotation period. They are accordingly readily excited by tidal interaction with a perturbing 
body when the characteristic time associated with the orbit is comparable to the rotation 
period, which is expected naturally when the rotation period and orbit become tidally cou- 
pled. They may then play an important role in governing the secular orbital evolution of the 
system. 

Inertial modes excited in close binary systems in circular orbit were considered by PP and 
Savonije & Papaloizou (1997). Wu (2005)a,b considered the excitation of inertial modes in 
Jupiter as a result of tidal interaction with a satellite and excitation as a result of a parabolic 
encounter of a planet or star with a central star was studied by Papaloizou & Ivanov (2005), 
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hereafter referred to as PI and Ivanov & Papaloizou (2007), hereafter referred to as IP. The 
latter work was applied to the problem of circularisation of extrasolar giant planets starting 
with high eccentricity. In that work the planet was assumed coreless. Ogilvie & Lin (2004) 
and Ogilvie (2009) have considered the case of a cored planet in circular orbit around a 
central star and found that inertial waves play an important role. 

The importance of the role played by inertial waves in the transfer of the rotational energy 
of a rotating neutron star to gravitational waves via the Chandrasekhar-Friedman-Schutz 
(CFS) instability was pointed out by Andersson (1998). Later studies mainly concentrated 
on physical mechanisms of dissipation of energy stored in these modes that limit amplitudes 
of the modes, and, consequently, the strength of the gravitational wave signal. In these 
studies either numerical methods or simple local estimates of properties of inertial modes 
were mainly used, see eg. Kokkotas (2008) for a recent review and references. 

An analytical treatment of problems related to inertial waves, such as eg. finding normal 
mode spectra and eigenfunctions, and coupling them to other physical fields, etc., is difficult 
due to a number of principal complicating technical issues. 

In particular, the dynamical equations governing the perturbations of a rotating body 
(called planet later on) are, in general, non-separable, for compressible fiuids. When such 
fluids are considered and rotation is assumed to be small, a low frequency anelastic approx- 
imation that filters out the high frequency modes is often used (see eg. PP). This simplifies 
the problem to finding solutions to leading order in the small parameter RlQ'^/{GM), where 
Q is the rotation frequency, G is the constant of gravity and M*, R^, are the mass and radius 
of the planet. In this approximation eigenfrequencies of inertial modes are proportional to 
Q, while the form of the spatial distribution of perturbed quantities does not depend on the 
rotation rate. However, even when this approximation is adopted, the problem is, in general, 
non-separable apart from models with a special form of density distribution, see Arras et al 
(2003), Wu (2005)a and below. 

Additionally, the problem of calculating the inertial mode spectrum and its response 
to tidal forcing is complicated by the fact that in the inviscid case the spectrum is either 
everywhere dense or continuous in any frequency interval it spans (Papaloizou & Pringle, 
1982). This is in contrast to the situation of, for example, high frequency p modes, which 
are discrete with well separated eigenvalues. When the anelastic approximation is adopted 
the singular ill posed nature of the inviscid eigenvalue problem is seen to come from the 
fact that the governing equation is hyperbolic and the nature of the spectrum is determined 
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by the properties of the characteristics (eg. Wood 1977). A discrete spectrum is beheved 
to occur when there are no such trajectories that define periodic attractors. Otherwise the 
inviscid spectrum is continuous. Then, when a small viscosity is introduced the spectrum 
becomes discrete but normal modes have energy focused onto wave attractors (see eg. Ogilvie 
& Lin 2004). Given these complexities it is desirable to work with and compare a variety of 
analytical and numerical approaches. 

Coreless inviscid rotating planets with an assumed spherical or ellipsoidal shape have 
a discrete but everywhere dense spectrum that makes difficulties for example with mode 
identification and application of standard perturbation theory. However, numerical work 
indicates that there are well defined global modes that can be identified and followed through 
a sequence of models (eg. Lockitch & Friedman, 1999, hereafter LF, and PI). In this paper 
we investigate the inertial mode spectrum of a uniformly rotating coreless barotropic planet 
or star and its tidal response by a WKBJ approach coupled with first order perturbation 
theory and compare its eigenvalue predictions with numerical results obtained by a variety of 
authors and find good agreement apart from some unidentified WKBJ modes that are near 
the limits of the spectrum and for which the perturbation theory appears not to work. For the 
identified modes we also find remarkably good agreement for the form of the eigenfunctions. 
This indicates that they can be represented at low resolution with small scale phenomena 
being unimportant, meaningful mode identification (in that the modes can be followed from 
one model to another) and at least first order perturbation theory works for these modes. 

This is also confirmed in a following paper (hereafter referred to as PIN) where we 
investigate the inertial mode spectrum and its tidal response by numerical solution of an 
initial value problem without the anelastic approximation. We are able to confirm the validity 
of the anelastic approximation and the applicability of the first order perturbation theory 
developed here for demonstrating this as well as estimating eigenvalues. Thus a suggestion 
of Goodman & Lackner (2009) that tidal interaction might be seriously overestimated by 
use of the anelastic approximation is not confirmed. 

A WKBJ approach to the same problem was also considered by Arras et al (2003) and 
Wu (2005) a. However, in this work only terms of leading order in an expansion in inverse 
powers of a large WKBJ parameter A (see the text below for its definition) were taken into 
account and treatment of perturbations near the surface and close to the rotational axis 
were oversimplified. As a consequence, although their results are correct in the formal limit 
A — > oo, they cannot be used to make a correspondence between WKBJ modes and those 
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obtained numerically, or an approximate description of modes with a scale that is not very 
small. In this paper we treat the problem in a more extended way, considering terms of the 
next 0(A^^) order together with an accurate treatment of perturbations near the surface 
and close to the rotation axis. Additionally, we consider a frequency correction of the next 
order, 0(A~^), for modes having non-zero azimuthal number, m. 

We checked results obtained with use of the WKBJ formalism against practically all 
numerical data existing in the literature finding good agreement in practically all cases. 
Therefore, we can assume that our formalism may be applied to provide an approximate 
analytic description of inertial modes, including those with large scale variations, where the 
WKBJ approach might be expected to be invalid. Also, different quantities associated with 
the modes may be described within the framework of our formalism or its natural extension, 
such as the tidal overlap integrals (see PI and IP), quantities determining the growth rate due 
to the CFS instability and decay of inertial waves due different processes, eg. by non-linear 
mode-mode interactions (see eg. Schenk et al 2002, Arras et al 2003). Thus, the formalism 
developed here may provide a basis for the analytic treatment of inertial waves in many 
different astrophysical applications. 

The plan of the paper is as follows. In section |2] we briefly review the basic equations and 
their linearised form for a uniformly rotating barotropic planet or star. In section [2731 we go 
on to consider these in the anelastic approximation which is appropriate when the rotation 
frequency of the star is very much less than the critical or break up rotation frequency. We 
give a simple physical argument why we expect this approximation to be valid in this limit 
even when the sound speed tends to a small value or possibly zero at the surface of the 
configuration. In section 12.51 we give a brief discussion about when discrete normal modes 
may be expected to occur such as in the case of a coreless slowly rotating planet with surface 
boundary assumed to be either spherical or ellipsoidal. We then present a formal first order 
perturbation theory that can be used to estimate corrections to eigenfrequencies occurring 
as either a consequence of terms neglected in the WKBJ approximation or the anelastic 
approximation. The latter application is tested by a direct comparison with the results of 
numerical simulations in PIN. Section [276] concludes with a brief account of the form of the 
anelastic equations in pseudo-spheroidal coordinates in which they become separable for 
density profiles of the form p oc (1 — r^/i?^)^, where r is the local radius, R^, is the surface 
radius and /3 is a constant. (Arras et al 2003, Wu 2005a). 

In section [3] we develop a WKBJ approximation for calculating the normal modes which 
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is based on the idea that in the short wavelength hniit these modes coincide with those 
appropriate to separable cases which include the homogeneous incompressible sphere as 
a well known example. Solutions of a general WKBJ form appropriate to the interior of 
the sphere are matched to solutions appropriate to the surface regions where they become 
separable which is the case when the density vanishes as a power of the distance to the 
boundary as is expected for a polytropic equation of state. This matching results in an 
expression for the eigenfrequencies given in section 13.51 

In section 13.61 we go on to develop expressions for the eigenfunctions appropriate to 
any location in the planet including the rotation axis and the critical latitude region where 
one of the inertial mode characteristics is tangential to the planet surface. These solutions 
are then used to obtain corrections to the eigenfrequencies resulting from density gradient 
terms neglected in the initial WKBJ approximation in section 13.91 In section |4] we compare 
the corrected eigenfrequencies obtained from the WKBJ approximation with those obtained 
numerically by several different authors who used differing numerical approaches and find 
good agreement even for global modes. A similar comparison with the results of numerical 
simulations for a polytropic model with positive results is reported in PIN. We also compare 
the forms of the eigenfunctions with those obtained in Ivanov & Papaloizou (2007) and find 
a good agreement even for global modes. 

Finally in section [STTj we discuss our results in the context of the evaluation of the overlap 
integrals that occur in evaluating the response to tidal forcing. We show that these vanish 
smoothly in the limit that the polytropic index tends to zero and we indicate that they 
vanish at the lowest WKBJ order and are thus expected to vanish rapidly as the order of 
the mode increases. We go on to summarize our conclusions in section 15.21 



2 BASIC DEFINITIONS AND EQUATIONS 

In this section we review the formalism and equations we adopt in this paper. As much of 
this has been presented in previous work (PI, IP) only a brief review is given here. 

In what follows we continue to investigate oscillations of a uniformly rotating fully convec- 
tive body referred hereafter to as a planet, focusing on the low frequency branch associated 
with inertial waves. 
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2.1 Framework for linear perturbation analysis 

The planet is characterised by its mass M*, radius /?* and the associated characteristic 
frequency 



where G is the gravitational constant. We adopt a cylindrical coordinate system (tz?, 0, z) 
and associated spherical coordinate system (r, 0, 9) with origin at the centre of mass of the 
planet. 

In this paper we make use of the Fourier transform of a general perturbation quantity, 
say Q, with respect to the azimuthal angle and the time t in the form 

Q = ^ ( exp(im0) / daQm exp{—iat) + cc j , (2) 

■Ul \ J —QC J 

where the sum is over m = and 2 and cc denotes the complex conjugate of the preceding 
quantity hereafter. The reality of Q implies that the Fourier transform, indicated by tilde 
satisfies Qml^") = Q*_„X~'^\ The inner products of two complex scalars Yi, Y^-, that are 
functions of ro, z and t are defined as 

{Yi\Y2)= f zudwdz{Y*Y2), (3) 

where * denotes the complex conjugate. Note that the definition of the inner product differs 
from what is given in IP where the planet's density p was used as a weight function. Integrals 
of this type are always taken over the section S of the unperturbed planet for which = 0. 



2.2 Linearised equations of motion governing the response to tidal 
perturbation 

We assume that the planet is rotating with uniform angular velocity fl. The hydrodynamic 
equations for the perturbed quantities take the simplest form in the rotating frame with z 
axis along the direction of rotation. 

Since the planet is fully convective, the entropy per unit of mass of the planetary gas 
remains approximately the same over the volume of the planet, and the pressure P can be 
considered as a function of density p only, thus P = P{p)- As the characteristic oscillation 
periods associated with inertial modes are in general significantly shorter than the global 
thermal timescale we may adopt the approximation that perturbations of the planet can 
be assumed to be adiabatic. Then the relation P = P{p) holds during perturbation as well 
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leading to a barotropic equation of state. In the barotropic approximation the hnearised 

Euler equations take the form (see PI) 

_| + 20x^ = -Vtr. (4) 

where 

W = cIp / p + ¥"' + ^'''\ (5) 

^ is the Lagrangian displacement vector, p is the density perturbation, Cg is the adiabatic 
sound speed, \l/™* is the stellar gravitational potential arising from the perturbations and 
\1/^^* is an external forcing potential, say, the tidal potential in the problem of excitation of 
inertial waves by tides, see PI and IP. 
The linearised continuity equation is 

p' = -V-(pO- (6) 

Note that the centrifugal term is absent in equation (jlj) being formally incorporated 
into the potential governing the static equilibrium of the unperturbed star. The convective 
derivative -^ = ^ as there is no unperturbed motion in the rotating frame. 

Although incorporation of the perturbation to the internal gravitational potential presents 
no principal difficulty, in this paper, for simplicity we neglect it, setting \&*"* = 0. This pro- 
cedure known as 'the Cowling approximation' can be formally justified in the case when 
perturbations of small spatial scale in the WKBJ limit are considered. However, it turns 
out that when low frequency inertial modes are considered the Cowling approximation has 
been found to lead to results which are in qualitative and quantitative agreement with those 
obtained numerically for global modes obtained with a proper treatment of perturbations 
to the gravitational potential (see below). Therefore, we do not expect that the use of the 
Cowling approximation can significantly infiuence our main conclusions. 

Provided that the expressions for the density and sound speed are specified for some 
unperturbed model of the planet, the set of equations (HI -TBI) is complete. Now we express 
the Lagrangian displacement vector and the density perturbation in terms of W with help 
of equations (jl]) and 1^, and substitute the result into the continuity equation from which 
we obtain an equation for its Fourier transform in the form 

a^AWrr, - aBWrn - CW^n = cr^d^i^^:' - W^), (7) 

where d = AQ"^ — cr^, and 
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B^J-^lL, c=-4tf^fp^). (9) 

to ow oz \ az j 

It is very important to note that the operators A, B and C are self-adjoint when the inner 
product 



{yV^\W2) = / dzrx7dwW^W2, (10) 

J y 

with V denoting the volume of the star. Here these operators are assumed to act on well 
behaved functions and the density is taken to vanish at the surface boundary. Also when 
m 7^ A and B are positive definite and C is non negative. When m = 0, A remains 
positive definite if consideration is restricted to the physically acceptable variations W that 
conserve mass, this constraint eliminating the possibility that W is constant. 

When the Cowling approximation is adopted equation ([7]) fully specifies solutions to the 
problem of forced linear perturbations of a rotating barotropic planet. In the general case, 
a complete set of equations is described in PI. 

2.3 The anelastic approximation 

When \^^* = 0, equation ([7]) leads to an eigenvalue problem describing the free oscillations 
of a rotating star in the form 

a^AWm - aBWm - CWm = -a^d^Wm, (11) 

Assuming that rotation of the planet is relatively slow such that the angular velocity 
fi ^ fi^,, these may be classified as / or p modes with eigenfrequencies such that \a\ > ^2* or 
inertial modes with eigenfrequencies |cr| ~ fi. The / and p modes exist in non rotating stars 
and can be treated in a framework of perturbation theory taking advantage of the small 
parameter a = \^/cr\ (see eg. Ivanov & Papaloizou 2004 and references therein). 

On the other hand for inertial waves a is of order unity, such a perturbation approach 
cannot be used. Since, in general, equation (ITT!) is rather complicated even for numerical 
solution, in order to make it more tractable the so-called 'anelastic approximation' has been 
frequently used (see eg. PP, Lockitch & Friedman 1999 and Dintrans & Ouyed 2001) for 
which the right hand side of ( ITTl) is neglected. 

In order to justify this approximation we note that for eigenfunctions that are non 
singular everywhere in the planet, we can crudely estimate the derivatives entering equation 
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(ITTl) as \VW\ ~ kW/R^ and |V^iy| ~ k'^W/Rl, where the parameter k > I. Consider first 
the interior region of the planet where we approximately have c^ ~ GM^jR^. It follows from 
equation (ITTll that the left hand side and the right hand side can be respectively estimated 
as ~ ^ pk^WmlR^^i and ^^pWml(?s- The ratio of these is of order 

~ < 1. (121 

This estimate is, however, not valid near the boundary of the planet where c^ — )■ and the 
left hand side of the inequality (1121) diverges. However, in the same limit the terms containing 
the density gradient on the left hand side of (TTT]) will dominate terms involving the second 
derivatives of Wn- Thus in this limit the magnitude of the contribution from terms on the 
left hand side of ([7]) may be estimated to be 

^' ^ kWjR, ~ VL^^lpkWmjcl, (13) 

dr 

where we remark that it follows from hydrostatic equilibrium that close to the surface 

\dpjdr\ ~ pGM^/{c^gRl). Accordingly, when r ^ R^ the ratio of the terms on the right 

and left hand sides of equation (fTTj) can be estimated as 

«1 (14) 



fi2A; 

From equations flT^ and f lT^ it follows that when Q '^Q^ the terms determining deviation 
from the anelastic approximation are small compared to the leading terms everywhere in the 
planet. Accordingly, in the slow rotation regime, we can use this approximation to find the 
leading order solutions for eigenfrequencies and eigenfunctions and then proceed to regard 
the terms on the right hand side of (ITTl) as a perturbation. 

The validity of the anelastic approximation in the context of the tidal excitation of inertial 
modes has been recently questioned in a recent paper by Goodman & Lackner (2009) on 
account of the divergence of the terms on the right hand side of equation lTTT!) as r — )■ i?^, 
although an actual demonstration of its failure was not given. In fact the above discussion, 
which also applies to equation (jT)) as this differs only by the addition of a forcing term, 
indicates that these terms are never important provided Q/Q^, is sufficiently small. This is 
to be expected because as Q is reduced, the structure of the modes remains unaffected in 
the anelastic approximation whereas the radial width of the region where terms on the right 
hand side of equation ( fTTT) might become comparable to any other terms shrinks to zero. 

We also note that the vanishing of the normal velocity at the boundary in the anelastic 
approximation is correct in the limit Q/Q^, — )■ as the ratio of the horizontal to normal 
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components there can be shown using the above arguments to also be on the order of 
{Q/Q^:)"^. Finally in PIN, we find by comparing the results of tidal forcing calculations using a 
spectral approach with the anelastic approximation, to those obtained using direct numerical 
solution of the initial value problem, that it gives good results even when Q/Q^: is not very 
small. 

2.4 Self-adjoint formalism 

It was shown by PI and IP that both quite generally and also when the anelastic approx- 
imation is used equations dTj) and (TTTj) can be brought to the standard form leading to an 
eigenvalue problem for a self-adjoint operator. Here we describe the approach, which leads 
to the self-adjoint formulation of the problem in the anelastic approximation. 

The self-adjoint and non negative character of the operators A, B and C is made use 
of to formally introduce their square roots, eg. A^'^, defined by condition A = A^"A^'^, 
etc. As is standard, the requirement of non negativity, makes the definitions of these square 
roots unique. The positive definiteness of A (see above discussion) also allows definition of 
the inverse of A^/^, A^^^^. 

Let us consider a new generalised two dimensional vector Z with components such that 
Z = {Zi, Z2) and the straightforward generalisation of the inner product given by equation 
f irUj) . It is now easy to see that equation ([7]) is equivalent to 

oZ = HZ^ S, (15) 

where 

_ f A-^/^BA-^'^ A-^'^C^'^ \ 
^ ~ \ C'/^A-^'^ j ' ^^^^ 

and the vector S has the components 

(A-i/Vrf4(^;^^*-w^j, 0). (17) 

Note that as follows from (llSp the relation between the components of Z and Wm can be 
taken to be given by 

Zi = aA^/^W^, Z2 = C^'^W^. (18) 



Since the off diagonal elements in the matrix flT6|) are adjoint of each other and the diagonal 
elements are self adjoint, it is clear that the operator H is self-adjoint. Equation ( 1T5|) can 
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be formally solved using the spectral decomposition of H. We now make a few remarks 
concerning the spectrum. 



2.5 The oscillation spectrum of a rotating fluid contained within an 
axisymmetric domain 

It has been known for many years (see eg. Greenspan 1968, Stewartson & Rickard 1969) 
that the eigenvalue problem we consider is not well posed in the inertial mode range —2Vt < 
a < 2VL. This is because in this spectral range the eigenvalue equation (fTTj) becomes a 
hyperbolic partial differential equation with boundary conditions specified on the planet 
boundary. The form of the spectrum depends on the behaviour of the characteristics, which 
correspond to localised inertial waves, under successive reflections from the boundary. Note 
that these reflections maintain a constant angle with the rotation axis rather than the normal 
to the boundary. The situation was conveniently summarised by Wood (1977) (see also Fokin 
1994a,b and references therein). There are three types of behaviour of the characteristic paths 
for frequencies in the inertial mode range. They may all close forming periodic trajectories, 
they may be ergodic, or there may be a finite number of periodic trajectories that form 
attractors. The first two types of behaviour are believed to be associated with discrete 
normal modes while the third type leads to wave attractors and a continuous spectrum. 
The homogeneous sphere within a spherical or ellipsoidal boundary exhibits the first two 
kinds of behaviour and has discrete normal modes which form a dense spectrum (eg. Bryan 
1889) while the same system with a solid core has wave attractors (eg. Ogilvie & Lin 2004, 
Ogilvie 2009). Note the characteristics behave in the same way for all spheres or ellipsoids 
with a continuous density distribution so that these should have normal modes. Note too 
that in the limit of very short wavelength modes only the second derivative terms matter 
in equation (TTT]) and the system becomes equivalent to the two dimensional case studied 
by Ralston (1973) and Schaeffer (1975). In that case the normal modes are associated with 
the frequencies for which all characteristic paths are periodic. They form a dense spectrum 
and are infinitely degenerate. From this discussion we expect the modes of a system with 
a continuous density distribution to approach the same form as those of the homogeneous 
sphere, an aspect upon which we build our later WKBJ approach. 
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2. 5. 1 Formal solution of ( [73]) in the anelastic approximation 

From the above discussion we expect the normal modes for the cases of interest to form a 
discrete but dense spectrum. The anelastic approximation can be implemented by setting 
Wm = in equation flT7|) . In this case we can look for a solution to flTSj) in the form 



Z = Y,akZk, (19) 

k 

— * 

where Zk are the real eigenfunctions of H satisfying 

(^kZk = HZk, (20) 

the associated necessarily real eigenfrequencies being ak- 
Substituting (120]) into ( 1T^ we obtain 

«fc = ^ ^1^ ^. 7- (21) 

The operator H induces the inner product and associated orthogonality relation for eigen- 
functions according to the rule 

< Zk\Zi >= akCri{Wk\AWi) + {Wk\CWi) = S^i, (22) 

where 

Wk = C-'/^Z^ = al^A-^l''Z\. (23) 

Using (TT7]) and (12T]) we explicitly obtain 

k Nk{a-ak) cf 

where 

Nk = al{Wk\AWk) + (WklCWk) (25) 

is the norm. The decomposition (II 9p should be valid for any vector F with components 
(F, 0), where F is any function of the spatial coordinates. The second component of this 
equality shows that in order for this to be valid an identity 

T.^^J^Z'2=0 (26) 

k ^^k 

must be hold (IP). This identity allows us to represent the relation ( 124]) in a different form 
(PI): 

Wm = <^dj: .r.f_^A Wk\-^,^TWk. (27) 

, l\k{0 (Jk) Cs 
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Note that in response problems such as the problem of excitation of the inertial waves during 
the periastron flyby, in order to take account of causality issues correctly when extending to 
the complex a plane, one should add a small imaginary part in the resonance denominator 
in (1271) according to the Landau prescription: (cr — a^) — )■ [a + iv — a^.), where z/ > is a 
small real quantity. 

2.5.2 Corrections to the anelastic approximation 

When external forces are absent and the potential ^^* is set to zero, equation ([7]) (or, alter- 
natively, equation ( !T5|l ) defines the full eigenvalue problem. Under very general assumptions 
it was shown by IP that this problem can be formally solved in an analogous manner. How- 
ever, it is rather difficult to use the general expressions obtained by IP without making 
further approximations. Here we note that, given that the spectrum is discrete, we may 
find conditions satisfied by the eigenf unctions and eigenvalues by replacing \&^* by —Wm in 
equations (12^ and (127|) . These conditions relate any eigenfunction, now equated to Wm and 
its associated eigenvalue a to the eigenf unctions and eigenvalues of the anelastic problem. 
Proceeding in this way we go on to form the quantity 

2 7 

aai{Wi\AWm) + {Wi\CWm) = - / ''\ {Wi\^Wm). (28) 

[a - ai) cl 

where Wi is an anelastic eigenfunction and we have made use of the orthogonality relation 
(El. 



As argued in section 12. 3[ the quantity on the right hand side can be regarded as a 
perturbation where the small parameter is e = (fi/f]*)^. Provided an eigenfunction can be 
identified as e — )■ and is non degenerate with W^ — ^ W^h it follows from fl28|) that in this 
limit 






i^-^i) = -^iWi\^^Wi), (29) 



where do = 4fi^ — af. 

The spectrum of inertial modes is dense. This may lead to a potential difficulties in 
identifying and following modes as parameters change as we discussed above. However, it 
is possible to argue that this problem can be alleviated for large scale global modes by for 
example modifying the eigenvalue problem by adding terms that have a very small effect on 
the global modes but spectrally separate close by short scale modes. Dintrans and Ouyed 
(2001) adopt such a procedure by adding a viscosity and this enables them to identify and 
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follow global modes. Note that a similar situation would result if conservative high order 
derivative terms were added that preserved the self-adjoint form of the problem. Numerical 
work presented below and in PIN also confirms that global modes have a clear identity and 
can be followed as parameters change provided that the angular frequency is sufficiently 
small. Thus we both expect and verify the validity of the expression (1291) in this limit. 
For larger values of Vt one should take into account a possibility of mixing between two 
neighbouring large scale global modes to explain results of numerical calculations, see PIN. 
In this case expression ( l29l) should be modified in an appropriate way. 



2. 5. 3 Eigenvalues corresponding to opposite signs of m 

In the next section we find solutions of the eigenvalue problem in the WKBJ approximation. 
It will be shown that the corresponding eigenvalues and eigenmodes are independent of the 
sign of m to two leading orders. This is explained by the fact that to that order solutions 
are determined only by operators containing second and first derivatives in equation ((Tj). On 
the other hand it follows from the same equation that the only dependence on sign of m is 
determined by the operator B which does not contain any derivatives of W . 

In order to find the first correction to the WKBJ eigenfrequencies that depends on sign of 
771, we treat the operator S as a perturbation. This leads to a change in the eigenfrequency 
that can be found by using the same formalism that lead to equation (l29ll but then simply 
replacing aidipWi/c^g in that equation by —BWi. Equation (!29|) then gives 

a-ai = -2mnaf I f dwdz-^wA /Ni. (30) 

Note that since g;^ < it follows from equation ( l30ll that when fi > the sign of a — o"/ is 
proportional to the sign of m. 

2.6 A form of equation ([7]) valid for a spherical planet 

In what follows we assume that an object experiencing tidal interactions can be approximated 
as having a spherically symmetric structure. In this case it is appropriate to use another form 
of ([7]) with \I'^* = 0, which is especially convenient for an analysis of WKBJ solutions. We 
can obtain this from (jTj) using the fact that for a spherical star -^p = ^-^p and -^p = f ^p. 
We obtain 



OZ'^ 



"^-Th 



' o dW , dW 

ozu oz 



2maVt-d{a/VtK{r)f W], (31) 
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where we set, for simplicity, W = Wm, A is the Laplace operator, H = — ^^ is a character 



istic density scale height and ^^(r) = J ^^ -^ , where M{r) is the mass enclosed within a 
radius r. Note that we use the hydrostatic balance equation —cj.^ = ^^ to obtain equa- 
tion ( 13T]) from equation ([7]). The last term in the second square braces on the right hand 
side describes correction to the anelastic approximation. It is discarded when the WKBJ 
approximation is used. 



2.6.1 Pseudo-spheroidal coordinates 

When the density approaches a constant value, H tends to infinity and the right hand side 
of equation ( 13T|) vanishes. In this case it describes an incompressible fiuid, see eg. Greenspan 
(1968). It was shown by Bryan (1889) that in this case this equation is separable in special 
'pseudo-spheroidal' orthogonal coordinates defined by the relations 



G7 = -R* 



{1-XJ){1-XI) R^XiX2 , ,, . , , <y .ooN 

z = , where the constant parameter /i = — — . [S2) 



\ (1 — ^^) fj, 2U 

Since the governing equations are invariant to the mapping (a, m) — )■ (— a, — m), without 
loss of generality we assume from now on that /i > for all modes while m can have either 
sign. Also, from equation (|3T|) it follows that the modes should be either even or odd with 
respect to the refiection in the equatorial plane z — )■ —z. Therefore, it is sufficient to consider 



only the upper hemisphere {z > 0, y/w^ + z^ < i?*). In this region we can assume that the 
variables xi and X2 are contained within the intervals [//, 1] and [0, /x], respectively. A detailed 
description of this coordinate system can be found in eg. Arras et al (2003), Wu (2005)a. 
Using the new variables equation ( l3T1) takes the form 

{D2 - D,)W = ^,^^ _^^,^^^ {Aw - {xl - xDim^i - V(l - /x2)(fi/n^(r))2)iy) (33) 



where 



d , r,d m 



2 



A = a;i(l - xl){xl - fi^)- ^2(1 - xl){xl - /i^)^— , (35) 

OXi 0x2 

and the quantities r, H, Q:^{r) are understood to be functions of the variables Xi and X2- 

It is easy to see that the eigenfunctions of the operators Di are the associated Legendre 

functions, P^{xi), and we have 

D,P:r{x,) = A2p-(x,), (36) 
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where A = Jh'{h' + 1). Let us stress that as the domains of Xi and X2 are not [—1, 1], z/ is 
not necessarily an integer. 

In some important cases equation fl33l) is separable. Firstly, when the gas is incompress- 
ible, the right hand side of (133|) is zero. In this case it follows from equation (!36|) that the 
solution can be represented as product of two associated Legendre functions. 

Secondly, as was mentioned by Arras et al (2003) and later explored in detail by Wu 
(2005) a,b when equation (!33|) is considered in the anelastic approximation it is separable for 
planetary models with density profiles of the form 

where C is a constant. These models include the incompressible one which corresponds to 
/3 = 0. It was also noted by Arras (2003) and Wu (2005)a,b that for polytropic models with 
equation of state 

p = kp^, (38) 

close to surface the density distribution has the form 

p = Dx", (39) 

where x = 1 — r/R^, n = 1/(7 — 1) is the polytropic index, and D is a constant. In the 
asymptotic limit r ^^ R^ this expression coincides with what is obtained from equation 
(!37|) with P = n. This proves that when polytropic models are considered equation (!33|) 
is separable in a plane parallel approximation often adopted close to the surface. Here we 
would like to note that this is valid even when the anelastic approximation is relaxed since 
close to the surface we have fix(^) ~ fi* = const and the additional term appearing in this 
case in the braces on the right hand side of (15^ has the same spatial structure of a term 
already present in the anelastic approximation. 

3 WKBJ SOLUTIONS FOR THE NORMAL MODES 

In general equation fl33l) should be solved numerically. We can, however, look for analytical 
solutions to ( l33l) in the WKBJ approximation assuming that solutions are fast oscillating 
functions in the planet's interior. The first and second derivatives of these functions are 
assumed to be proportional to first and second power of a large parameter A, the value of 
which is specified below. This problem has been analysed before by Arras (2003) and Wu 
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(2005)a who obtained expressions for eigenvalues and eigenfunctions for the problem of free 
oscillations in the inertial mode spectral range. Here we revisit the problem, taking into 
account terms that appear at the next order in an asymptotic expansion of the quantities 
of interest in a power series in 1/A. This will allow us to obtain analytic expressions which 
agree with numerical results, even for the rather small values of A, appropriate to global 
modes (see below). 

3.1 Natural units 

In what follows in order to simplify notation we express all dimensional quantities in natural 
units. These are such that the spatial coordinates, density, angular velocity and sound speed 



are expressed in units of i?^,, the mean density p = 3M^:/{AtcRI),JGM^:/RI and JGM^,/R^ 
respectively. All other quantities of interest are expressed in terms of powers of these basic 
units. 

3.2 WKBJ solutions 

It is easy to find from either f l3T]) or fl33|l that in the planet's interior far from the rotational 
axis, the WKBJ solution should have the form 

WwKBj = -^(Cie*^'^+("+) + C2e^^<^-("-) + cc) + 0{\), (40) 

^/pW A 

where 0±(m±) are arbitrary functions of 

u± = w± z, (41) 

the constancy of which defines the characteristics of equation fl33l) . Acceptable forms for the 
functions 0± have to be determined by matching the solution (140|) to approximate solutions 
valid near the surface boundary and near the rotational axis. It turns out that this matching 
is possible if the WKBJ solution has the form 

WwKBj = ^^(cos(A|/i + 0i) cos(Ay2 + 02)), (42) 

where yi^2 = arccos(xi^2) and 0i 2 are constants to be determined. One can readily check with 
help of the coordinate transformations (l32l) that this form agrees with the general expression 
( HOj) . see equation ( l85l) below. Since yi^2 are multivalued functions of xi,2 we should specify 
a one-valued branch of these. Taking into account that our calculations will be done for 
positive values of xi,2, we assume below that values of 1/1,2 are in the range [0, vr/2]. 
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For simplicity, in the main text we are going to consider the modes even with respect to 
reflection z — )• —z, called hereafter 'the even modes'. For example, such modes are excited 
by tidal interactions since tidal potential is an even function of z. The case of the modes 
odd with respect to this reflection ('the odd modes') can be dealt with in a similar way. This 
case is considered in Appendix A. 

From equation (1321) it follows that reflection of the coordinate z leads to the reflection 
of the coordinate X2 such that X2 — )■ —X21 while the coordinate 2/2 changes according to the 
rule 1/2 — ^ TT — y2- We readily find that P2l) is unchanged under this transformation provided 
the phase 

02 = -7rA/2, (43) 

(see also eg. Wu 2005a). We remark that the same result is obtained by requiring that the 
derivative of fH2l) with respect to 1/2 vanish on the equator where ^2 = 7r/2. 

3.3 Matching near the rotation axis 

In the WKBJ approximation sufficiently far from the rotational axis all terms proportional 
to W give small corrections to the solution (I40p and are formally discarded. However, when 
tu — )■ and, accordingly, Xi — ?■ 1, it follows from equation (1M|) that the term proportional 
to rn^W in the expression for the operator Di diverges in this limit and should be retained. 
When this is done the phase 0i can be found from condition of regularity of W close to the 
rotation axis ro = 0. 

We begin by using the WKBJ solution already found to develop an approximate expres- 
sion for W that is appropriate for small values oi 5 = 1 — Xi and which can be matched 
at large distances from the rotation axis. An appropriate expression for W which can be 
matched to the correct WKBJ limit sufficiently far from the rotational axis is 

W oc -^^^^^-pjCosA(y2 - vr/2)iy„(5), (44) 

where we take into account that the factor y/w entering (HOj) is proportional to the product 
{I — x'lY/'^il — xlY^"^ oc 5^/'^ {I — x'lY/'^, see equation fl32|) . and the factor 5^^^ is formally 
incorporated in the definition of Wa{5) which is to be found by imposing the condition of 
regularity on the rotation axis. In order to do this we obtain an equation for Wa from equation 
fl3T|) (or fl33l) ) that retains terms containing the derivatives and terms that potentially diverge 
in the limit cc — )■ while other terms can be discarded. 
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From equations (13T]) and (!33l) it follows that Wa{5) satisfies equation (!36|l in the limit of 
small 5 



The solution to fH5|) regular at the point 6 = can be expressed in terms of the Bessel 
function 

W,^J\„^\{V2\6'/'), (46) 

where we assume from now on that A is positive!. In the limit of large {X6) the asymptotic 
form of the expression ( H6|) is 

Wa oc (5-1/^ COS (V2\6'/^ - 1^1 1 - l) • (47) 



It is easy to see that when 6 is small yi ~ V^S. Therefore, from equations ( l44l) and ( l47l) it 
follows that the solution has the required form (H2]) provided that 

0i = -|"^l2-4' (48) 

and we have, accordingly, 

1 „ „ /, , ,vr 



W 



Wi^'BJ 



^ ^F, F = cos (^\yi - |m| I - ^ cos (^Xy2 - A0 . (49) 

Note that the phase f H5]) . which can in fact be verified with reference to the incompressible 
sphere, differs from that given in Arras et al (2003) and Wu (2005)a. This disagreement 
is due to an oversimplified treatment of the WKBJ solution close to the rotational axis in 
these papers. 

3.4 Matching at the planet surface 

The eigenvalues appropriate to the problem of free oscillations can be found by matching 
the solution f H9|) to approximate solutions valid near the surface of the planet. In pseudo- 
spheroidal coordinates ( 132|) the equation determining the upper hemispherical surface of the 
planet {r = R^,z > 0) has two branches: 1) Xi = fi,0 < X2 < fi and 2) fi < Xi < 1,X2 = fi- 
In order to simultaneously consider solutions to equation ( l33l) that can be close to either of 
these branches, we introduce two new coordinates 6i, with i = 1 corresponding to the first 
branch and i = 2 corresponding to the second branch, that are defined by the relation 

Xi = fi±6i, (50) 

^ In the approximation we consider the final expressions are independent of the change of sign of A. 
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where the sign + (— ) corresponds to the 1st (2nd) branch, and assume later on that the Si 
are small. 

The form of the solutions close to the surface depends on the density profile. In what 
follows we the consider the planet models with a pol5d;ropic equation of state for which the 
density profile close to the surface is given by equation ( l39l) . The variable x = 1 — r/R^, 
entering equation ( l39l) can be expressed through Xj and 6i as 



..^M_4^, (51) 

where we assume from now on that the upper (lower) sign corresponds to the 1st (2nd) 
branch and the index i takes on the values 1 ( first branch) and 2 (2nd branch) with j ^ i. 

We now look for solutions close to the surface that have X6j large but X6i small. This 
is possible because in the WKBJ theory A is a large parameter. The domain for which 
X6i is small for both i = 1 and i = 2 is called the critical latitude domain and will be 
considered separately below. Solutions valid in all of these domains must match correctly 
on to a solution of the form (B9l) in order to produce a valid eigenfunction. 

Using equations (13 9 p and fH2]) we can look for a solution close to the surface in the form 

W oc \x^j - /i2|-"/2(i _ a;2)~i/4 cos(^^. + (f).)Wi{5i), (52) 

where we also use equation (132|) in order to express the factor tu"^/^ in terms Xj setting 
6i = there. Substituting this expression in equation ( l33l) and taking the limit (Jj — )■ we 
obtain 

d^ T.r d ,,, 1 

, W + n — W -\ 

d6^ d6 (l-/i2 

where, for simplicity, we omit the index (i) in the quantities 6i and Wt, and 

m, = m-4^(l-/i2) — . (54) 

We recall that the term proportional to Q/Q^, gives the correction to the anelastic approx- 
imation. Since in the low frequency limit fi/fi^, is assumed to be much smaller than unity, 
this term is small and we approximately have m^, ^ m. 

Equation ( l53l) can be brought into a standard form by the change of variables 

C = 2iKd, ^ = e'^^W = e^/^iy, where k = , (55) 

VI -/U^ 

Adopting these we obtain 

C,^" + (n - 0^' -\{n^ %X)^ = 0, (56) 



'^TH^ + ^37^ + 7^ ^(A25±nm,)l^ = 0, (53) 
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where a prime denotes differentiation with respect to C, and x = n'm^/{\^l — fi'^). This is the 
conffuent hyper-geometric equation. Its solution that is regular at the surface is expressed 
in terms of the confluent hyper-geometric function $(a, b, z) as 

Vroce-^/^$((n±ix)/2,n,C). (57) 

Note that this solution is similar to solutions of the Schrodinger equation with the Coulomb 
potential describing wave functions belonging to continuous part of its spectrum, ( see eg. 
Landau & Lifshitz 1977). 

In the limit of |fi;5| ^ 1 we obtain from (1571) 

'' - <-"'' (fw^W) '^^ (' (-^ ^ t '■' 1^1 + t) ) + -) • (««' 

where G{z) is the gamma function. Since the quantity x is assumed to be small we can 
approximately write 

l/n^in T ^X)) ^ (1 ± #(|)f )/r(|) ^ exp (#(|)|) /r(|), (59) 

where il}{z) = -^T{z) is the psi function. In the same approximation equation (158|1 can be 
rewritten in the form 

W oc C^"/^ cos (Kl ± I (^in Id - ^(^)) - ^) . (60) 

After substituting the result expressed by equation (!60|) into (1521) the resulting expression 
should be of the general' form (given by H2|) evaluated close to the surface. This, however, 
cannot be realised on account of the presence of the factor In |d in (160!) . This term, having 
a coordinate dependence of order of 1/A after removing a constant phase would formally 
require terms of that order that are not accounted for in the expressions (HOll and (H2ll to 
enable matching, therefore to the order we are currently working, it is discarded. Since only 
this term depends on the sign of m and on the correction to the anelastic approximation, 
both dependencies are absent in the resulting approximation. 

Another way of obtaining solutions to (153|) compatible with the form (142|) inside the 
planet is to set to zero the small quantity x oc 1/A in equation fl57|) . In this case the solution 
can be expressed in terms of a Bessel function such that a 

iy(5)oc5(i-")/Vi/2(„-i)(«:<5). (61) 

^ In order to obtain equation H61|l we use the relations ^{a, b, z) = e^<I>(6 — a, b, —z) and ^{a, 2a, 2z) oc z^'^~'^e^ Ja-i/2{i^)^ 
see eg. Gradshteyn & Ryzhik 2000, pp 1013, 1014. 
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Note that this expression is equivalent to (!57|) when the anelastic approximation is adopted 
and m = 0. When K(5 — )■ oo we get 

iy,oc5r/^cos(-^=i=5.-^), (62) 

where we the index (i) has been restored and we use the exphcit expression for k. Substituting 
fp2]) into equation (15^ . taking into account that the factor (|a;^ — /i^|(5i)~"/^ oc y/p, and that 
close to the surface we have 

Hi ^ arccos/i =F ,^ ' ., ■ (63) 

3.5 Determination of the eigenfrequencies 



It can now be seen that the expression ( 1521) has the required form (l42l) provided that the 
phases (pi satisfy appropriate appropriate conditions. However, these phases have already 
been determined from the requirements of regularity on the rotation axis and symmetry 
with respect to reflection in the equatorial plane and are accordingly specified through 
equation (l49ll which equation(!52l) must match. 

It is readily seen that the expressions (H9l) and ( l52l) can be compatible only for particular 
choices of A and /i. These compatibility conditions determine the eigenspectrum of the 
problem in the WKBJ approximation. They are easily found from equations fj49|) . f l62|) and 
( l63|) to be given by 

—n + nki = A arccos(/i) |?Ti| and —n + 7rA;2 = -^ ( arccos(/i) ) . (64) 

Here fci and ^2 are positive or negative integers that must be chosen in a way which ensures 
that the angle arccos(/i) belongs to the branch < arccos(/i) < |. 
Adding the above relations we obtain 

A = 2/ + n+ |m| + -, (65) 

where I = ki + k2. Substituting ( l65l) into the first expression in ( l64l) we obtain an expression 
for the eigenfrequency 

a ( fc+H/2 + (ri + l)/4\ 
/. = - = cos [. ^ j , (66) 

where we set k = ki from now on. 

As shown in Appendix A the modes with different symmetry with respect to reflection 
z —7- —z (both the 'even' and the 'odd' modes) can be described by the same expression fj66|) 
provided that the expression for A changes to 
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\ = p + n + \m\ H — , (67) 

where the integer p is even for the even modes while for those with odd symmetry p is odd. 
For the WKBJ approximation to be vahd A should be large, and, accordingly, 1^1. 
We would like, however, to consider all values of / and k allowed by our assumption that A 
is positive and arccos(/i) belongs to the interval (0,7r/2). These conditions imply that / is 
positive and lead to inequality: 



[I + n) \m\ 
+ 



< k< 



n 
^+4 



(68) 



4 2 

where [Q] means that integer part of Q is taken. 

When / and k are sufficiently large one may neglect other quantities in the argument of 
the cosine in equation (l66ll . In this case one gets /i = cos(7rA;/(2/)) - an expression obtained 
in previous papers (see Arras 2003 and Wu 2005a). One may also consider the limit of an 
incompressible fluid by setting n = in (l66l) . In this case the expression ( l66l) gives the 
correct asymptotic eigenfrequencies appropriate to the high order modes of pulsation of an 
incompressible fluid in a rotating spherical container, see Appendix B for details. 

3.6 A general expression for eigenfunctions close to the surface of the planet 

The purpose here is to establish an expression for W that is approximately valid in the whole 
region close to the surface where the separation of variables is possible and the eigenfunction 
can be written as the product of functions of Xi alone and X2 alone. Also, this expression 
should approach the WKBJ expression (H9!) in the limit of sufficiently large Ax in order to 
have the norm that will be given by equation ( 1911) below. 

Close to the planet surface we have x = 1 — r ^ 1, the density profile can be represented 
in the form (!39l) . and equation ([7]) becomes separable in the pseudo-spheroidal coordinates 



(xi, X2) (Arras et al 2003, Wu 2005a). As described already in section [3^ in these coordinates 
the region close to the surface is described by two branches Xj = /i ± 5j, [i = 1,2) see 
equation ( l50l) . We denote these branches as the (+) branch for which i = 1 and the (— ) 
branch for which i = 2, respectively. 

When one of the coordinates, Xj, say, is sufficiently far from the value Xj = fi the eigen- 
function is proportional to the expression given by equation ( l52|l . In practice the requirement 
that Xj is far from /i is that \XSj\ be large. When this parameter is of order unity or less, 
Xj is considered to be close to /i. When both coordinates are close to /i in this sense, the 
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eigenfunction is proportional, with, in the hmit of large A, proportionality factor being slowly 
varying, to the product Wi{6i)W2{S2) , where Wi{xi) can be found from equation ( 16T1) . From 
equation (l32ll it follows that when xi ~ X2 ~ /i the spherical polar angle 6 is close to the 
critical latitude defined by 6' = arccos(/i). Accordingly, we describe this region as the region 
near the critical latitude. 



3. 6. 1 An expression for the eigenfunction near the critical latitude 

In order to obtain an expression for the eigenfunction that is valid near the critical latitude we 
proceed as follows. At first we consider a region of the planet sufficiently far from rotational 
axis, where A(l — Xi) ^ 1. We start from the form of solution given by equation fl52l) in the 
form 

W oc (±(/i2 - x^ i)-"/2(l - x2 i)-V4 cos(A7/2,l + </'2,l)W^(5l,2), (69) 

where here and below the first index and upper sign (second index and lower sign) correspond 
to the (+) branch ((— ) branch). The function W{y) satisfies equation (l53ll with ??i* set to 
zero. As we discussed above the term proportional to m^, gives a correction which will 
be calculated below. The fact that the function W{y) is normalised in order to have the 
appropriate limit in the case of large y is stressed by the overbar. We have 

W{y) = y^y(^-")/V(„_i)/2(«:2/), (70) 

where we recall that k, = A/-\/l — fi'^. 

Equation fl69|) is not valid near the critical latitude where both xi and X2 are close to /x. 
To obtain a modified form that is valid, the function cos(A|/2,i + 02. i) must be replaced by a 
function that matches this when A|?/2,i ~" arccos(/i)| is large but which takes has the correct 
form to result in the proportionality of the eigenfunction, as indicated above to W{6i)W{62), 
when this is small. It can be seen that an expression for W having the required properties 
can be written in the form 

W± = D;'/\1 - xl,)-'/\±{fi' - xl,)^'/'W{6^,2)A^2!lW{A2,l), where (71) 



A2,i = ±Jl - /u2(|/2,i - arccos(/i)). (72) 
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Here the quantity D^ = -D((l — /i^)/i)~"', where D is defined in equation ( l39l) B We remark 
that in the hmit k,6i^2 — ^ oo together with the hmit k,Ai^2 — j- oo we can use the asymptotic 
expansion of Bessel functions 



where 

, ^ (-1)^ r(f + 2fc) ^ (-1)^ r(f + 2fc + i) 

'^ 22'=(2A;)!r(f-2A;)' '^ 2(2'=+i)(2A; + 1)! r(f - 2A; - 1) ' ^ ^ 

together with equation fl64l) to show that expression ( 17T|) attains the required form specified 
in equations (|42]) and (iH]) . 



3.6.2 An expression for the eigenfunction near the pole and surface 

In the region close to the pole of the planet we have 0:2 ~ jU and 6 = 1 — a^i ^ 1. The 
discussion given above excluded consideration of this domain which needs to be considered 
separately. Close to the pole but away from the surface the solution is given by equation 
(I44p ■ From very similar considerations to those above, close to the surface where \X62\ is 
of order unity or less and to the pole where 1 — xi ^ 1, the solution is proportional, to 
within, in the limit of large A, a slowly varying proportionality factor, to the product of 
W{S2) and the solution given by equation fH6|) which is valid near the rotation axis. Thus 
in this domain W oc Wa{S)W{62). An an approximate solution for Wpoie that satisfies the 
required matching conditions and also reduces to the form (17T1) in the limit A(l — xi) 3> 1 
can be written as 



W^poze = {-irD;'/\i - xl)-'/\xl - ^^r-/'^j^JUXy,)W{S2), (75) 

where the factor (—1)*''^ takes into account that when ki is odd Wa and Wwkbj differ by 
sign in the matching region, see equation f l64l) . 



3.7 An approximate expression valid over the whole surface domain 

We now use an interpolation procedure to combine expressions derived above, that are valid 
in separate domains inside the planet, to form single expressions that can be used over 
the whole domain. To do this we introduce a function 77 [x] that is defined for x G [0, 1] 
and belongs to C*^°°\ We stipulate that 77 [x] decreases monotonically with x and is such 

^ Let us stress that the A2.1 is not necessarily small contrary to 52,1 defined through equation II50I I. 
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that ri[x\ = 1 when x G [0, 1 — x*] and ri[x] = when x G [x^,, 1] The quantity x^, is a 
parameter such that 1/2 < x* < 1. An exphcit form of r][x] as well as the value of x=k are 
not important for our purposes. This is because all representations in the different domains 
attain a matching asymptotic form in the planet interior and elsewhere. 

3. 7. 1 An expression for W- 

We may now write down an approximate solution valid close to the surface for all Xi G [/i, 1] 
for W- which we denote as W^ as 

Wl = Ti[z^]W^ + (1 - r][zi])Wpoie, (76) 

where zi = {xi — /i)/(l — /x) and we recall that W^ in the above is to be obtained from 
equation f lTTj) . The expression fl76|) can be rewritten in another useful form, which explicitly 
shows that the solution is separable close to the surface: 

Wl = D;'/\1 - xl)-'/\xl - /i')-"/'iyi(xi)iy(52), (77) 

where 



m = v[^i]K^'W{A,) + (1 - ^[^i])^!^j|„|(AyO. (78) 

3.7.2 An expression for W+ 

In this case we formulate an expression valid close to the surface and for all X2 G (0,/i). 
Close to the equatorial plane and away from the critical latitude, and K,{fi — X2) = KA2 ^ 1 
it is convenient to represent the function W^ in terms of an asymptotic series in ascending 
powers of (/tA2)~^. Substituting the series ( j73l) in ( ITTI) for W+ we obtain: 

W,, = D;'/\1 - xl)-''\^? - xl)-^I^W,,{x,)W{5,), (79) 

where 

_ 00 4 00 TD 

W,,{x,) = (-1)'=^ cos Xiy, - -) E ^-^ + (-1)(^^+^) sin A(y2 - vr/2) ^ ^^^-^y (80) 

where the coefficients A}, and Bj. are given in equation ( FM|) and we take into account that 
fi;A2 — rni /A = X{y2 — 7r/2) + 7ik2, according to equations f l64|) and f l72|) . 

On the other hand it is convenient to use the expression (ITTI) directly in the region 
close to the critical latitude. The expressions ( ITT]) and ( I79|) can be combined with help of 
the function r][x] to provide a single function W^, which is analogous to the function Wl 
discussed above and can be used for X2 G (0,/i): 
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Wl = D;'/\1 - xl)-'/\^^' - xlr-/'W,{x2)W{6^), (81) 

where 

W2{X2) = v[z2]A^2^'W{A2) + (1 - v[z2]mg{x2), (82) 

and Z2 = I — X2//X. 

3.8 Calculation of the mode norm 

In order to calculate different quantities related to a particular mode we need an expression 
for the mode norm given by equation (125|) . One can show that when A is large enough the 
interior of the planet gives the dominant contribution to the integrals determining the norm. 
Thus the expression (H9|) for the eigenfunction is appropriate. In addition one can simplify 
the expression ( l25l) by noting that eigenf unctions satisfy equation (I7j) with the right hand 
side set to zero. Furthermore, the term involving the operator B in equation ([7]) may be 
neglected. This is because it does not contain second derivatives and therefore contributes 
a higher order correction to the WKBJ approximation as discussed above. Additionally, for 
the same reason, only the terms proportional to the second derivatives in the operators A 
and C need to be retained since these terms give the leading contributions to the norm 
being proportional to A^. Thus we have 

a'^AWwKBj ~ CWwKBj, (83) 

and 

N ^ 2{WwKBj\CWwKBj) = Sn^ 1^ dwdz (^f\ , (84) 

where in order to evaluate the norm ( 125|) , we use the explicit form of C given by equation ([9]) 
and, after an integration by parts, adopt equation (H9l) for the eigenfunction. The quantity 
F can be expressed in the form 

F = ^(cos(Ay+ - VI/+) + cos(Ay_ - ^^)), (85) 

where 

y± = y2±yi, ^± = ^|m| + ^ ± ^A. (86) 

Note that y±. can be readily expressed in terms of coordinates (xi, X2) and (cc, z): 



cosy± = a;ia;2 =F V 1 — ^iv'^ ~ x^ = fiz ^ (wl — ij?)w, (87) 
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where we use equations (132|) . From equation (1871) it follows that the quantities y± are constant 

on characteristics of equation l^^, see equation ( 14T1) . By differentiating ( l85l) we then obtain 

dF ^ _fiX / sin(Ay+ - ^+) ^ sin(Ay--^-) \ 
92; 2 Y sin|/_|_ sin|/_ y 

In the limit of large A the integral of square of the derivative fj88|) in f l84|) can be approx- 
imately calculated taking into account that only average values of sin^(A|/± — \l/±) ~ 1/2 
give a significant contribution to the integral. In this way we obtain 

where 

/ = / rdrd9{ — . — + — . — } = 2%. (90) 

Jv sin y+ sin ?/_ 

Note that in order to evaluate the integral ( l90l) we use the fact that from ( 1871) we have 

sin^ y± = 1 — r^cos^(^ ± arccos(/i)). The integral is elementary and most easily done by 

noting that it is easily shown to be independent of /i, and can accordingly be evaluated 

setting /i = 1. Substituting fl90|) into (|8^ we obtain a very simple expression for the norm: 

N = 2n{jjXny. (91) 

3.9 Calculation of the frequency correction cr™ 

As follows from equation (!66|) the eigenfrequencies of the modes are degenerate with respect 
to change of sign of m in the approximation we use. The correction to the eigenfrequencies 
accounting for the dependence on sign of m, a™, is calculated above, see equation (l30l) . As 
follows from this equation, a^ is determined by the integral 

I = - f dwdz-^Wl (92) 

Jv dw ^ ^ ^ 

The integral fl92l) has contributions from the interior of the planet where the mode 
eigenfunction is given by ( H9|) and also close to the surface where equations (I7T]) and (!75|) 
apply. Contributions to the integral close to the surface arise from both the (+) and (— ) 
branches. Accordingly, we have / = /^^■^■^ + /+ + /^. At first let us evaluate the contribution 
from the inner region, j'^^^J _ In order to do this we substitute ( l49l) into ( l92l) . thus obtaining 

jWKBj ^_J drdOF^^ In p. (93) 

Since the quantity F is rapidly oscillating we use the average value of F^, < F"^ >= 1/4 in 
(I93D. Hence 
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A Jo dr "^ 4 p^ ^ ^ 

where pc is the value of the central density of the planet and p* ^ pc is a value of density 
close to the surface of the planet (r = r,, = 1 — e) above which the contribution to / has 
to be determined by a separate treatment of the surface region. With help of equation ( 139|) 
equation ( p3|) can be rewritten in the form 

/^™ = ^(lng-nlnx,), (95) 

where x* = 1 — r^ = e is the dimensionless distance from the surface corresponding to the 
density p*: p* = p(a;*). 

Now let us evaluate the contribution to integral from the region close to the surface. For 
definiteness, let us consider the (+) branch where 6i = Xi — p is assumed to be small. An 
analysis of contribution of the region close to the critical latitude to the integral fl92|) shows 
that this contribution is small and can be neglected. Thus, we can use an expression of the 
form f l52|) . choosing i = 1 and j = 2 there, that correctly matches ( H9|) in the interior. We 
substitute this in ( l92l) and change the integration variables from cylindrical coordinates to 
pseudo-spheroidal coordinates, using the fact that 

zudzudz = - — -^dxidx2 ~ -, :;^dSidx2, (96) 

(see eg. Wu 2005a). Also, we use the approximate density profile given by equation ( l39l) and 
express the radial variable x there in terms of 6i and X2 with help of ( 15T1) . In this way we 
obtain 

/+ =n j dx2d6i{l - xlY^/^Fid'l-^W^di), (97) 

where F2 = cos(A|/2 + 02)- The integral ( 1971) should be evaluated in a volume bounded by 
the surfaces a; = and x = x^,. 

First let us evaluate the integral J„ = J d5i5i~^W{5iY. Using the explicit expression 
flTO]) for W{y), we obtain 

In = l£*dyjf„_,y,iy). (98) 

Here we note that the integral converges only for the case n > 0. The integration variable 
y = K,6i with y^: = k,5^. The condition y = y^, defines the surface x = x^,. From equation flHT]) 
we have 

5. = ^^-/\^* . (99) 

p^ -xi 



/+ = -(ln(5„ft:(l - /i2)/^x,) / dx2{l - xlY^'^ - <l>+(/i)), (102) 

2 Jo 
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The integral (!98|) logarithmically diverges when y* — )■ oo. When y^ is sufficiently large it 
can be represented in the form 

J„ = ^ln(5„y,), (100) 

where the constant B^ can be calculated numerically for a general value of n|j • Substituting 
(llOOj) to (jnZl), using the averaged value of F2, < F^ >= 1/2, we obtain 

/+ = - / dx2{l - xlY^I'' ln(5„y,), (101) 

2 JO 

where a multiplicative factor of two has been applied in order to account for the contributions 
from both the upper and lower hemispheres. Using equation ([22]) we can bring (llOip to the 
form 

M 


where 

$+(^) = r dx{\ - x')-^/2 ln(/i2 _ x^). (103) 

h 

The same approach can be used to evaluate the integral corresponding to the (— ) branch 
with the result: 

I = -(ln(5„fi:(l - /i2)^x,) / dxiil - x?)-i/2 - $-(/!)), (104) 

2 J ^ 

where 

$-(^) = f dx{l-x^)-^^^\n{x'^-ij.^). (105) 

J fj, 

Thus, the integral corresponding the the region close to the surface, /^ = J"*" + I~ can be 
evaluated as 

I' = -(- ln(5„fi:(l - /i2)/ia;,) - <^tot), (106) 

where 

$^„^= / dx{l-x^)-^^^ln\x^- fx^\. (107) 

Jo 

It can be shown (eg. Prudnikov et al 1986) that the last integral does not depend on fi: 
^tot = — vrln2. Substituting this value to (11061) . remembering that the total integral I = 
jWKBj _|_ js ^^^ adding (I106p to (1951) . we obtain the final expression for the integral 



/ = ^ln 

4 



g(4i?„A/iv/l-/iT 



(108) 



* The factor ^ is determined by the form of asymptotic expression ofJ^(j/) fs ./ — cos (j/— ■|-J'— t)- Substituting this to l|98|l , 
integrating the result between two values of 3/ = 1/1,3/2, then assuming that j/2 3> 3/1 we obtain this factor. 
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Note that the integral f llOSp does not depend on position of the matching point x*. Substi- 
tuting it and the norm ( l9Ti) to the expression for the frequency correction ( 130|) we have 



< = — In 



Pi 
D 



(45„A/iv/l - /i2)- 



(109) 

In Appendix B we show that the expression (11091) has a correct form in the hmiting case of 
an incompressible fluid ra = 0. 

The expression (I109p is not valid when |/i| is sufficiently close to 1. We recall that one 
can prove that the absolute value of any eigenfrequency must be less than or equal to 20, ( 
see eg. PP). This condition may be violated when AyU = 1 — /i ^ 1 and the correction 
(11091) is added to the unperturbed frequency (!66l) . A similar constraint may be obtained 
from consideration of the assumptions leading to (1 1001) . Indeed, the matching radius should 
obviously be smaller than the radius of the star, and, accordingly, x^, < 1. For small values 
of A/i this condition together with equation (!99|) leads to 

5, < /i, (110) 

where we assume that a 'typical' value of X2 entering (1991) is of the order of ~ yU ~ 1. On 
the other hand, for the validity of the asymptotic expression (!97|l we should have K(5* > 1, 
and therefore 

S. < ^. (111) 

Combining inequalities (IllOp and (lllip we obtain 

A/i > §, (112) 

where a coefficient C > 1 can be obtained from a more accurate analysis. Thus, our simple 
approach is likely to be invalid for eigenfrequencies with absolute values sufficiently close 
to 2Q. Therefore, we discard unperturbed values ( l66l) that lead to eigenfrequencies with 
absolute values larger than 2Q when the correction (I109p is added. 

The same analysis can be used to calculate the quantity a^" determining the correction 
to the anelastic approximation. As follows from equation f l29|) the correction is proportional 
to an integral very similar in form to that given by equation (!92|) . This integral also has 
contributions from the surface and the interior where the standard WKBJ solutions may be 
used. These contributions are comparable and so they should be matched at some radius. 
In fact, it may be shown that the surface contribution is equal to the surface contribution 
to f l92|) given by equation (I106p . The internal WKBJ contribution is more complicated and 
should, in general, be evaluated numerically. Therefore, for simplicity, we adopt a different 
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Figure 1. Positions of the eigenfrequencies calculated in the WKBJ approximation and the numerical results obtained by LF 
on the a/Q, axis are plotted for the case \m\ = 2. The open circles and squares are obtained from the WKBJ approach and 
the pluses and stars give the results of the numerical computations of LF. The circles (pluses) correspond to modes having 
I = linin virhilc the Open squares (stars) correspond to modes with / = Imin + 1- 

approach when deahng with the correction to the anelastic approximation. We calculate 
0"°" numerically using the expression ( l29l) for several 'global' eigenmodes (i.e. modes with 
a large scale distribution of perturbed quantities). Such modes are of especial importance 
in applications of the formalism. For example, as discussed in PI, IP and in PIN, those 
mainly determine the dynamic tidal response of the planet. Since o"^" oc A^^, we expect that 
corrections to small scale WKBJ modes are less significant. They are, therefore, neglected. 



A COMPARISON OF ANALYTICAL AND NUMERICAL RESULTS 



In this section we compare the frequencies obtained from the approach described above with 
those obtained by a number of authors who have employed a variety numerical methods. 

To obtain the values from the above analysis we use equation fl66l) for the eigenfrequen- 
cies, and, in the case of non zero m, add the expression for the correction to the frequency 
given by f llOQp . We compare the results for polytropes with polytropic indices n = 1, and 
1.5. The quantities i?„ entering (11091) for these cases were obtained numerically with help 
of equation (llOOp . We found Bi ^ 14.4 and -B3/2 ~ 5.97. The range of allowed values of k 
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Figure 2. As for Fig. [Tlbut the left panel illustrates the case m = and the right hand side panel corresponds to |m| = 1.. 
In the left panel we also show the locations of the WKBJ modes corresponding to I = 4 (diamonds) and the locations of two 
eigenfrequencies calculated by Dintrans & Ouyed (2001) (triangles). 



O +0 O 04 



□ *□ m D 



□ * D 



, M M M M M M M M M M M M M M M M M M M 
-2 -1.8-1.6-1.4-1.2 -1 -OS-0.6-04-O2 02 04 06 0.8 1 1.2 1.4 1.6 1.8 2 

a/a 



O O O O O + O C) 



nn*n )C nae n*n[] 



-2 -1.8-1.6-14-1.2 -1 -08 -0.6 -0.4 -02 02 0.4 0.6 08 1 1.2 1.4 1.6 1.8 2 



Figure 3. As for Fig.[T]but for the case of |7Tt| = 3 (left panel) and \m\ = 4 (right panel). 



is found from equation fl68|) . Additionally, we shall discard modes, which have \fi\ close to 
or 1 (see discussion below). We comment here that it follows from equation (11091) that the 
frequency correction a^ becomes undefined for such values of /i. 
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Figure 4. The result of comparison of the odd modes. Circles, squares and diamonds show positions of the WKBJ modes, 
corresponding to 1° = 0, 1, 2, respectively. Pluses, stars and crosses give positions of frequencies numerically calculated by LF. 
The cases of |m| = 1 (left panel) and |m| = 2 (right panel) are shown. 
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Figure 5. As for Figs.|4]but for the case of \m\ = 3 (left panel) and \m\ = 4 (right panel). 

Most of the eigenvalues considered are for the model with n = 1, which has been exten- 
sively discussed, (see eg. LF and Dintrans & Ouyed, 2001). The amount of attention paid 
to this model is due to the fact that it approximately reproduces the density distribution 
of a cold coreless Jupiter mass sufficiently below the planet's surface. Some eigenvalues for 
planetary models (mostly the so-called 'global modes' with m = 2 and / = 1) were calculated 
by IP. The properties of the global modes in this case are very similar to those of a polytrope 
with n = 1. Since the different numerical approaches essentially agree with each other for 
the case n = 1 with non zero values of ni, we compare our results with the numerical values 
obtained by LF. 
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LF consider both even and odd modes and classify the mode order by an integer Iq, 
which is related to integer p defined in equation ( 1671) through Zq = p+\m\ — l. Let us remind 
that when the even modes are considered p = 21, see equation (l65l) . and for the odd modes 
p is related to the integer 1° classifying the odd modes as p = 2/° + 1, see equation (]A3|) of 
Appendix A. Using this definition, it is easy to see that LF give eigenfrequencies of modes 
having even symmetry with / = Z^m = 1 and / = 2 for the case |r7i| > and / = /^m = 2 
and / = 3 for the case m = 0. For the case of m = 1, we also considered the next order 
even modes with / = 4 and compared eigenfrequencies with what was obtained by Dintrans 
& Ouyed (2001). We also compare the WKBJ odd modes for m = 1,2,3 and 4 with the 
results of LS. The integer 1°, in this case, takes the values 1° = Z^j„ = 0, 1 and 2. For the case 
n = 1.5 the comparison of the WKBJ eigenfrequencies is made with the numerical results 
obtained by spectral methods in PI. 

The results of the comparisons for n = 1 are shown in Figs. [1]- [5], where positions 
of the eigenfrequencies within the allowed range —2 < a/Q < 2 are shown. The WKBJ 
eigenfrequencies, a are found from 

a = a, + <, (113) 

where a* = 2fi, and /i, a™ are given by equations fl66|) . (I109p . respectively. In our analytical 
investigation we assumed that the quantity fi, and, accordingly, a*, is positive, but allowed 
the sign of the azimuthal number m to be either positive or negative. For the purpose of 
this section it is convenient to take a different but equivalent point of view and assume that 
the sign of m is fixed: m > 1 and allow the quantities a* to have either sign. The positive 
and negative signs correspond to prograde and retrograde mode propagation with respect 
to direction of rotation of the planet, respectively. 

Let us first discuss the result of comparison of the even modes, which is presented in 
Figs. [TB 

In Fig. [1] illustrates the case with m = 2. This is the most important value when one 
considers the problem of dynamical tides excited by a perturbing companion. The WKBJ 
eigenfrequencies are shown as open circles for / = 1 and squares for / = 2, respectively. The 
corresponding numerical values are indicated by crosses and stars. As seen in Fig. [U the 
numerical and analytical values show quite good agreement with each other. This agreement 
is good even for the modes with the smallest possible value of / = 1, which have a global 
distribution of perturbed quantities over the volume of the star and therefore might not be 
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expected to be in any kind of agreement with the results of a WKBJ theory. This might be 
explained by the rather large value of the parameter A for these modes, being equal to 5.5. 
From the discussion above, its inverse A~^ ~ 0.18 is assumed to be a small parameter in our 
WKBJ expansions. There are, however, three unidentified modes for the case / = 1 as well 
as for 1 = 2. These modes have frequencies concentrating near the borders of the allowed 
frequency range with a ~ ±2 as well as in the region close to a = 0. It is possible that 
although these correspond to global modes for the incompressible {n = 0) model, these do 
not retain their character when n is increased to 1 and beyond through a strong coupling to 
short wavelength modes nearby in eigenfrequency (see discussion in section |23|) . Discarding 
these modes, hereafter referred to as unidentified modes, the relative difference between the 
analytical and numerical results is of order of or smaller than 10 per cent. 

In Fig. [2] we illustrate the case with m = in the left panel and the case with m = 1 in 
the right panel. For m = the eigenfrequency distribution is symmetric with respect to the 
reflection a — )■ —a, and, therefore, only positive values of a are shown. As for the previous 
case, when m = the agreement between the WKBJ values classified as 'physical' and the 
results of the numerical study is quite good, even when the 'global' mode with a ~ 1.3 is 
considered, the relative difference being of the order of 3 per cent for this mode . Again, 
this may be accounted for by the relatively large value of A = 5.5 for ?7i = and Imin = 2. 
Additionally, on this plot we show the results of a calculation by Dintrans & Ouyed (2001). 
They calculated five modes for m = 0. Three of them have locations nearly the same as 
those obtained by LF. They are, therefore, not shown in the plot. As seen from the plot, the 
other two modes may be identified with with the analytical modes having I = 4. Note that 
the agreement gets better with increasing WKBJ order , /, as expected. 

When |?Ti| = 1 the agreement is similar to the case with |r7i| = 2 apart from the mode 
with / = 2 and k = 2, where the WKBJ value a ~ 0.33 is approximately twice as small as the 
numerical one. The reason for this disagreement is unclear to us. Note, however, that there 
seems to be a disagreement between numerical methods in this case. Only one of the modes 
with / = 2 and |m| = 1 obtained by LF can be reliably identified with an eigenfrequency 
given by Dintrans & Ouyed (2001). However, the 'global' modes with \m\ = 1 and / = 1 as 
well as all LF modes obtained for the case m = have their counterparts in the results of 
Dintrans & Ouyed 2001. 

In Fig. [3] we show the cases of relatively large values of \m\ = 3 ( left panel) and \m\ = 4 
(right panel). The agreement gets worse with increasing \m\, and in the case of the 'global' 
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mode with |r7i| = 4, / = 1 and the numerical value a/Vt ^ 0.6, the disagreement is of the 
order of 30 per cent. This may be explained by the fact that in our theoretical scheme the 
value of |m| is assumed to be much smaller than the value of A. However, for the mode with 
the largest disagreement the ratio |r7i|/A ~ 0.53 is not very small. 

In Figs, m and [5] the comparison of the odd modes is made. The results are similar to 
the previous case. Apart from the presence of unidentified WKBJ modes, all numerically 
obtained frequencies have well identified WKBJ counterparts. The agreement is getting 
better with increase of the WKBJ order 1° and is getting worse with increase of the azimuthal 
number m. Note a rather good agreement between the 'global' WKBJ and numerical modes 
corresponding to 1° = 0. In fact, as was shown by PP (see also Papaloizou & Pringle 1977), 
eigenfrequencies and eigenfunctions of these modes can be calculated analytically giving very 
simple results 

a = -^^, W = zw^. (114) 

m + 1 ^ ^ 

Note that neither eigenfrequencies nor eigenfunctions depend on the planet's density distri- 
bution in this case. 

In summary, we point out that agreement between the numerical and WKBJ frequencies 
is unexpectedly good taking into account the fact that the WKBJ theory should not, strictly 
speaking, be applied to modes with such small values of A. Apart from the existence of the 
unidentified modes in the WKBJ scheme and the two 'physical' even modes and the global 
odd mode corresponding to m = 4 with the rather large disagreements alluded to above the 
agreement between analytical and numerical results is of the order of or smaller than 15 — 20 
per cent for all the remaining identified 27 even and 35 odd modes. As we shall see below, 
good agreement is also found when the spatial distribution of the modes is compared. 

4.1 Properties of eigenfunctions 

In order to calculate distributions of the quantity W over the volume of the planet we 
use equation ( H9|) in the bulk of the planet and equation (I7m77p close to the surface and 
smoothly interpolate between the two regions. Since the function r][x] defined in section [3171 
is inconvenient for a numerical implementation we consider instead of it a function, which 
is zero in the regions [0, 1 — x*] and [x^,, 1] and represented as a ratio of two polynomials 
of X in the intermediate region, which are chosen in such a way to ensure that several first 
derivatives are equal to zero in both points x = 1 — x* and x = x*. 
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Figure 6. Distribution of W over the surface of the star for the n = 1.5 polytrope. The upper plots are obtained by numerical 
methods and the lower plots from the WKBJ theory. The values of eigenfrequencies are given in the text. From left to right 
the integers I and k determining the WKBJ eigenfrequencies are / = 1, fc = 0, i = 1, fc = 1 and I = 2,k = 2 respectively. 

Let US stress that for self-consistency we use the frequency a* (or /i) as given by equation 
fl66|) in those equations even when the frequency correction is non zero. As above, the 
numerical results for the n = 1.5 polytrope are taken from PI. The WKBJ results for n = 1 
are compared with those obtained for a realistic model of a planet of one Jupiter mass. This 
model has a first order phase transition between metallic and molecular hydrogen, which 
has been discussed, eg., in IP. 



A comparison of the different results is shown in Figs. [6] and [3 Note that in all cases 
shown in this section the same contour levels are used for the numerical and analytical data. 

In Fig. Owe show the distribution of W over the planet's volume for the m = 2 modes and 
n = 1.5. Numerical results are presented in the upper plots, which are taken from PI. These 
are for modes with a = —1.06 ( upper left plot), a = 0.67 (upper middle plot) and a = 0.67 
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Figure 7. As for Fig. \E\ but for m = 0. The upper plots represent numerical results for a planet with a realistic equation of 
state, the lower plots are calculated from the WKBJ theory assuming that n = 1. The eigenfrequencies are given in the text. 



(upper right plot). The modes with a = —1.06 and a = 0.67 are the so-called two main global 
modes, according to PI. They mainly determine transfer of energy and angular momentum 
through dynamic tides induced by a parabolic encounter. The respective WKBJ counterparts 
have the smallest possible WKBJ order / = 1. The corresponding analytic eigenfrequencies 
are a = —0.99 (a* = —1.32) and a = 0.64 (cr* = 0.39). The distribution shown on the upper 
right plot may be identified with a next order mode having a = 0.435 and a^, = 0.29. One can 
see that there is a surprisingly good agreement between the analytical and numerical results. 
In particular, the retrograde mode represented on the left hand side plots has a 'spot' in 
distribution at the angle ~ 7r/4 with respect to the rotational axis. This agrees with position 
of the critical latitude since arccos/i = arccos |(T*/2| ^ 0.277r. The distributions shown on 
the middle and right plots correspond to prograde modes. They have a well pronounced 
approximately vertical isolines. The main global mode may be distinguished from the mode 
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Figure 8. Distributions of W for the global odd modes are shown. The azimuthal number m = 1, 2, 3 from left to right. The 
upper plots represent the analytical results given by equation l|114|l . the lower plots show the corresponding WKBJ counterparts. 



corresponding to the next order by the number of nodes in the horizontal direction, this 
being one in the case of the global mode and two for the next order mode. We have checked 
that similar agreement exists between the WKBJ and numerical results corresponding to 
n = 1. Since the distributions are quite similar they are not shown here. 

For m = we compare the WKBJ results with calculations done by a spectral method 
for a model of a planet of Jupiter size and mass in Fig. [3 As in the previous case the 
upper plots correspond to the numerical results. From left to right the numerical values 
of the eigenfrequencies are 1.35 (the main global mode), 1.01 and 0.79. Their analytical 
counterparts have I = 2,k = 1, a = 1.33; I = 3,k = 2, a = 1.03 and I = 4,k = 3, a = 0.83 
respectively. Note that a more pronounced disagreement in eigenfrequencies corresponding 
to the mode represented on the right hand side plot is mainly determined by the fact that 
this mode has a distribution concentrated near the surface of the planet, where the equation 
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of state differs from that of a n = 1 polytrope. One can see that again there is very good 
agreement between the results. This is especially good for the main global mode represented 
in the plots on the left hand side. The agreement gets somewhat worse moving from right to 
left. This may be explained by a number of factors such as inaccuracies of the numerical and 
analytical methods as well as the physical effects determined by changes of the equation of 
state in the outer layers of the planet and the presence of the phase transition. These factors 
mainly influence modes with a small spatial structure while the large scale main global mode 
is hardly affected by them. 

Finally we consider the global odd modes and compare the analytic distributions given 
by equation flll4p with the corresponding WKBJ distributions for m = 1,2 and 3 in Fig. 
IHl Although there is a disagreement in position of the spot close to the critical latitude, 
which is situated on the planet's surface in the case of the exact analytic solutions and 
slightly interior to the surface of the planet in the case of the WKBJ distributions, there 
is a similarity in the distributions in the planet's interior. This is quite surprising since in 
this case the analytic distributions do not depend on the planet's structure at all while the 
WKBJ distributions are determined by the density distribution close to the planet's surface. 



5 DISCUSSION 

5.1 Overlap integrals 

As we pointed out in the Introduction, integrals of the form 



Qk = ^j%- with Q, = (i^ty,|$), (115) 



''- with Qfc = ( — 

k ^s 

where Wk corresponds to a particular eigenmode and $ is some smooth function, appear 
in astrophysical applications of the theory developed in this paper. In particular, as was 
discussed in PI and IP, integrals of this type enter in expressions for the transfer of energy 
and angular momentum transferred during the periastron passage of a massive perturber. 
These apply to the case when the spectrum of normal modes is discrete and they involve 
integrals of form (11151) . where $ = P^r"^ with P^ being the associated Legendre function. 
Assuming that Wk varies on a small spatial scale while the function $ is smoothly varying 
such integrals may, in principal, be evaluated using our formalism with help of a theory of 
asymptotic evaluation of multidimensional integrals, see eg. Fedoryuk (1987), Wong (1989). 
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However, some important integrals of form flllSp require an extension of our formalism, 
which can provide a smooth matching of the solution close to the surface to the WKBJ 
solution in the inner part of the planet that is valid at the next orders in inverse powers of A. 
This is due to cancellations of leading terms in corresponding asymptotic series. Since this 
problem appears to be a rather generic one we would like to discuss it here in more detail 
for the important case when $ = P2r'^ = Stu^. The overlap integral of this type determines 
excitation of the m = 2 modes which are the most important for the tidal problem (eg. PI, 
IP and see also PIN). Explicitly, we have in this case 

Qk = 3ldv(to^-^Wk), (116) 

where dV = dzvjdxu. Note that this integral must converge to zero in the incompressible 
limit n —7- as in this case it is well known that inertial modes are not excited in the anelastic 
approximation. This fact, however, is not obvious for the integral written in the form (I116P 
because close to the surface we have 

4~nCx("-^\ (117) 

with the constant C converging to a nonzero value as n — )■ 0. Therefore, as the eigenf unctions 
are regular, the integral Qk/n has a logarithmic divergence at the surface of the planet as 
n — i- 0. This raises the possibility that the overlap integral might converge to a nonzero value 
or behave pathologically as the incompressible limit is approached. 

In order to show that, in fact, this is not so and Qfc(n — )■ 0) — ?■ in a smooth manner, 
let us consider some fiducial models having the property that the quantity 

-I = -^^ (118) 

rp dr 

is constant. For models in hydrostatic equilibrium under their own gravity, constancy of uJq 
implies that ratio M{r)/r^, where M{r) is the mass interior to the radius r, is constant. The 
model must accordingly be incompressible. Goodman & Lackner (2009) obtained a wider 
class of models in hydrostatic equilibrium under a fixed quadratic gravitational potential. 
Because the potential is fixed independently of the mass distribution and there are no con- 
straints on the equation of state, such models may be constructed for an arbitrary density 
distribution. 

Now let use consider the integral 

Qfud = ?>jdV (w^^ulwA . (119) 
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For the fiducial models described above this is identical to the overlap integral (I116P where we 
note that we may adopt natural units such that the constant Wq, which should be identified 
with the surface value of GM{r)/r^ is equal to unity in that case. More generally the 
integrand in f lllQp can be transformed using equation of hydrostatic equilibrium f lllSp such 
that 

Qfud = -^jdV (^^wi\ =-3jdV (^§^Wk) . (120) 

Now let us consider equation flTTj) for free normal modes in the anelastic approximation 
by setting a = (Xk and the right hand side of this equation to zero. Then we multiply it by w"^, 
set m = 2, and integrate over dV. After removing derivatives of Wk by integrating by parts, 
assuming that the density vanishes at the surface boundary, it is easy to see that it follows 
from ( ITTl) that Qfud = in the anelastic approximation provided ak 7^ 2Q. This means, in 
particular, that m = 2 inertial waves cannot be excited in the Goodman & Lackner (2009) 
models in this approximation as was found by these authors when compressibility was fully 
taken into account (see PIN for additional discussion). 

Using the fact that Qfud = we may rewrite (11161) for models under their own self-gravity 
quite generally, adopting natural units, as 

Qk = Qk- Qfud = 3 1 dVw^^ M - ^^ j Wk. (121) 

Taking into account that the factor in the brackets is proportional to x for small x, and, 
accordingly, in this limit the integrand is proportional to nx^ we see that now the logarithmic 
divergence of Qk/n disappears and, therefore, it is clear from the representation (I12ip that 
the overlap integral indeed smoothly tends to zero in the limit n — )■ 0. 

The theory of asymptotic evaluation of integrals of the form (I12ip tells that the values of 
such integrals are determined either by inner stationary points, where gradient of the WKB J 
phase vanishes or contributions close to the surface or other parts of the integration domain, 
where the WKBJ approximation is not valid. From the expression of the Wk in the WKBJ 
regime fHOl) it follows that there are no stationary points in the inner region of the planet. 
Considering the regions close to the surface it appears to be reasonable to assume that the 
leading contribution is determined by the region close to the critical latitude, where a 'hot 
spot' is observed in distributions of Wk, see the previous section. In this region the quantities 
5i = Xi— fi and 62 = fi — X2 are small. We can use them as new integration variables in (I12ip 
with help of fl96l) . separate the contribution of this region to the integral by introduction of 
the functions ?7[(5i,2] defined in section 3.7 in the integrand, and decompose the quantities in 
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front of Wk in powers of 81^2 taking into account that x oc 8182 and w"^ ^ (1 — /x^) in the 
leading order. Assuming that Wk ~ W {8i)W {82) , where W{ii) is given by equation ( 170l) . it 
is easy to see that the leading contribution to f lllSp is given by a symmetric combination of 
two integrals involving Bessel functions 

Qk (X h{8i)l2{82) + l2{8i)h{82), (122) 

where 

h{y) = f t/yr/[|/]y"J(„_i)/2(«:y), Uv) = I rfyr/[y]|/("+i)j(„_i)/2(«:i/), (123) 

Jy=0 Jy=0 

where it is assumed that ri[y > y^] = 0, y* lies within the range of integration and l/n ^ 
y^, <^ 1. As was shown by Larichev (1973) the integral Ii{y) = for any particular form 
of the function ri[y]. Thus, the leading order contribution to the overlap integral from the 
surface region close to the critical latitude is equal to zero. In principal, one can look for 
the next order terms. However, in this case our simple approach to the problem seems to 
be inadequate since eg. the assumption that Wk can be represented as a product of two 
functions separately depending on the coordinates may be broken at this level, etc.. A more 
accurate approach is left for a future work. We note, however, that this cancellation means 
that the overlap integral should decay rapidly with increasing A, possibly being inversely 
proportional to a large power of A. This may qualitatively explain why a small number of 
relatively large scale modes are significantly excited by dynamic tides, see PI, IP and PIN. 

5.2 Conclusions 

In this paper we have developed a WKBJ approximation, together with a formal first order 
perturbation approach for calculating the normal modes of a uniformly rotating coreless 
planet under the assumption of a spherically symmetric structure. Matching of the general 
WKBJ form valid in the interior to separable solutions valid near the surface resulted in 
expressions for eigenf unctions that were valid at any location within the planet together 
with an expression for the associated eigenfrequencies given in section 13.51 Corrections as a 
result of density gradient terms neglected in the initial WKBJ approach were also obtained 
from formal first order perturbation theory. 

The corrected WKBJ eigenfrequencies obtained using the WKBJ eigenf unctions were 
compared with results obtained numerically by several different authors and found to be in 
good agreement, away from the limits of the inertial mode spectrum where identifications 
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could be made, even for modes with a global structure. We also compared the spatial forms 
of the eigenfunctions with those obtained using the spectral method described in IP finding 
similar good agreement. 

This is consistent with the idea that these global modes can be identified and that first 
order perturbation theory works even though they are embedded in a dense spectrum. 

In further support of this, the formal first order perturbation theory developed here is 
subsequently used to estimate corrections to the eigenfrequencies as a consequence of the 
anelastic approximation and is then compared with simulation results for a polytropic model 
with n = Im PIN. These different approaches are found to be in agreement for small enough 
rotation frequencies and also indicates that, as implied by the simplified discussion in section 
12.31 of this paper, that corrections as a result of the anelastic approximation are never very 
significant for the models adopted. 

Our results show that the problem of finding eigenfrequencies and eigenvalues of inertial 
modes allows for an approximate analytical treatment, even in the case of modes having a 
large scale distribution of perturbed quantities. 

Although we consider only the case of a polytropic planet, our formalism can be applied to 
a much wider context. Indeed, the approach developed here is mainly determined by the form 
of the density close to planet's surface, where we assume that it is proportional to a power of 
distance from the surface. Thus, we expect that our main results remain unchanged for any 
density distribution, which is approximately power-law close to the surface. In particular, 
according to our results, all models of type having approximately the same behaviour of the 
density distribution close to the surface should have approximately the same eigenspectrum. 

The formalism developed here can be extended for an approximate analytic evaluation of 
different quantities associated with inertial modes, such as overlap integrals characterising 
interaction of inertial waves with different physical fields, growth rates due to the CFS 
instability and decay rates due to various viscous interactions and non-linear mode-mode 
interaction. It may provide a basis for a perturbative analytic analysis of more complicated 
models, such as realistic models of star and planets flattened by rotation or models of 
relativistic stars. 

As we discussed above, for a given value of WKBJ order, Z, some modes are identified 
with modes obtained numerically while others remain unidentified. Eigenfrequencies of the 
unidentified modes are always either situated close to the boundaries of the frequency range 
allowed for inertial modes, a = ±2Q or situated close to the origin a = 0. We believe that 



Inertial waves in rotating bodies: a WKBJ formalism for inertial modes 47 

our theory is not applicable to these modes, and they develop a small scale contribution 
controlled by the closeness of the positions of their eigenfrequencies to a = and ±2^2, 
and thus effectively move to higher order than allowed for. Accordingly we de not consider 
these modes when comparing our results with results of direct numerical calculations of the 
excitation of inertial waves due to a tidal encounter reported in PIN. 
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APPENDIX A: THE EIGENFREQUENCIES OF THE 'ODD' WKBJ 
MODES 

In order to find the eigenfrequencies of modes odd with respect to reflection z — )■ —z the 
phase 02 should be chosen in such a way that the corresponding eigenfunctions are equal to 
zero at the planet equatorial plane, and, accordingly, W{z = 0) = W{ii2 = f ) = 0. From 
this condition and equation fH2l) we obtain 
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02 = -A| + ^ + TTS, (Al) 

where s is an integer. Analogously to the case of even modes the phase 0i is determined 
by equation (HHI) and the form of solution close to the surface is determined by equations 
( l62l) and ( l63l) . From these equations and equation ( lAlj) we get the compatibility conditions 
analogous to conditions 



— n + TTKi = A arccos u \m\ and — ra + ttKo = A arccos u tts, A2 

4 ^'^^ 2' ' 4 4 V2 ^ V 2 ^ ^ 

and, solving (lA2p for A and /i we find that 

A = 2r + |m|+n + -, (A3) 



where 1° = ki + k2 + s, and that the expression for fi is given by equation fl66l) . where (IA3p 
should be used. Comparing equations ( 165|1 and (jASp we see that both even and odd modes 
can be described by the same expression for A provided that 

X=p+\m\+n + -, (A4) 

where p = 21 for the even modes and p = 21" + 1 for the odd ones, respectively. 



APPENDIX B: THE OSCILLATION SPECTRUM OF AN 
INCOMPRESSIBLE FLUID CONTAINED IN A ROTATING SPHERICAL 
CONTAINER IN THE WKBJ LIMIT 

In this appendix we show that the eigenvalues obtained from our WKBJ analysis agree with 
the corresponding values for an incompressible fluid contained in a rigid spherical container 
in the WKBJ limit. 

Bl Eigenvalues for the incompressible case in the WKBJ limit 

It is well known that the spectrum of normal modes for an incompressible fluid in a ro- 
tating spherical container can be found analytically, ( see eg. Greenspan 1968, p. 64). The 
eigenfrequencies are determined from the equation 

-mPH(/.) = (l-;.^)^i^, (Bl) 

where P^{fi) is a Legendre function, s is an integer, we recall that /i = cr/(2fi) and note 
that although it is inconsequential for the use of flBJl). Greenspan's definition of m has the 
opposite sign to that used in this paper. 
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We determine the spectrum in the WKBJ hmit from (1B1|) . In this hmit s is a large 
parameter. For our purposes it is important to retain all terms of order of 0{s~^) and 
larger. We use the asymptotic form of Legendre functions in the limit of large s in the form 

Pf (cos(/)) oc -^=f== i cos ( ( s + - ) - - + m- ) + (m^ - - ) \f ^ z\ J^ a. \ -(^2) 



v/sm0 [ VV 2/"^ 4 2/ V 1) 2(s + f)sin 

Here we have used the well known properties of gamma functions to transform the asymptotic 

expression given in Gradshteyn & Ryzhik (2000) to the form ( 1B2I) . 

We now substitute f lB2p into (]B1|) after setting \x = cos0 in that equation. We then 

discard all terms < 0{s^^), thus obtaining 

f f 1\ . </* I i</'\ ■ 1 ■ f f 1\ . ^ I i^r 

— mcos s H — \- \Tn\— ^ ssm0sm [ s -\ — \- \m\ — 

\\ 2/^ 4 ' '2/ ^ VV 2/^ 4 ' '2 

1 // 1\ TT 7r\ 1 2 1 // 3 TT TT 

+ -cos(^(^s--j0-- + |m|-j + -(m - - sm((s + -)0 + - + |m|-). (B3) 



The quantity multiplying s in the first term on the right hand side of (IB3P should be 

close to zero in order that this first term be of the same order as the other terms on the 

right hand side. This condition gives 

, + l/4-W/2 + ^ 
^ s + 1/2 ' ^ ^ 

where j is an integer and 5 is a small higher order correction. Setting it to zero we have 

/ j + 1/4- \m\/2\ , , 



This expression agrees with equation (1661) when we consider the limit of incompressible fluid. 
To do this we set n = in ( l65l) and ( l66l) . assume that j — \m\ = k and set 

A = s + i. (B6) 

Comparing (]B6P with ( 1651) we obtain 

s = 2/+|m|. (B7) 

One can show that the requirement that s — |m| is an even number determines eigenmodes 
for which W has even symmetry with respect to the reflection z — ^ — 2;, see eg. Greenspan 
1968, p. 65. Substituting (IB4p in (IB3P we get an expression for the correction 5 



8 = _!^i±i(i±£H. (B8) 

Now we substitute f IBSP into f lB4l) to obtain 

/^ = /io + ^(^ + m^) + ^. (B9) 
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The last term 

77? 

f^T = Y, (BIO) 

gives the leading order difference in eigenvalues belonging to eigenf unctions with values of 
m of opposite sign. It is shown below that the expression (IBIOI) agrees with (11091) provided 
that (11091) is evaluated in the limit 77 — )■ 0. 

B2 The limit 77 —t- of the expression ( 11091) for the frequency correction a^ 

If 77 is directly set to zero, the integral fl98|) diverges at 7/ = 0. In order to find a limiting 
expression for f ll09p it is necessary to consider n to be very small and take the limit n — ;■ 
in the final expression only after potentially divergent terms have been cancelled. Thus we 
assume that 77 > but very small and integrate ( 19 8 p by parts to obtain 

^n = o / dyJl{y) = -y.J^iy*) -n dyy-J^—Ju, (Bll) 

2 JO 2 JO dy 

where u = {n — l)/2 and we remark that for 77 > {yJ^{y)) — > when ?/ — )■ 0. We use the 
well known relation 

y-rMy) = ^Ju - yJu+i, (B12) 

dy 

to eliminate the derivative of the Bessel function and obtain 

ry* 
nin = TT[y*Jui.y*) + dyyJy{y)J^+i{y)\. (B13) 

JO 

It is easy to see that the expression on the right hand side does not diverge when 77 — )■ 0. 
However, the integral entering the right hand side of (1B13P diverges at large values of y^. In 
order to deal with this divergence we use the asymptotic expression for the Bessel function, 
J^{y), valid at large values of its argument in the form, 

.,,)., ,„./X{eo.(.-.^-9-i-I™.„(.-.|-^)}. ,BM, 

From equation (lB14p it follows that the boundary term in (]B13p can be evaluated as 

y*Juiy*) = ^ (1 + cos {2y, - ^i)) + O^y*'^- (B15) 

Here we remark that because y* ~ Xx^,, where x* is the dimensionless distance to the surface 
and A is large, y* may be large when x* is small corresponding to being close to the surface. 
Thus use of the asymptotic expansions of Bessel functions for large values of their arguments 
can be justified. In addition, setting u = {n — l)/2 and z/ + 1 = 77 + 1/2, flB14p gives 

JLi{y)J^n±i{y) ~ — (sin (2y - 77^) + f] , (B16) 

2 2 7iy \ \ 2 / 2y I 
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and 



fy* 1 / 7v\ n 

I*iy*)= dyyJl^{y)J'Ui{y) — cos 2y, - n- + — Iny,. (B17) 

JO 2 2 ZTT V I J ITX 

Now, restoring complete precision, we can represent the integral entering (IB13P as 

y* 

dyyJu{y)Ju+i{y) = h{y*) + C„, (B18) 



where the quantity C„ does not contain any divergences and can be evaluated in the limit 
ra — )■ and y^ — > oo. 

Setting ra = in (IB17P and using the exact expressions for the Bessel functions given by 



J_.=^Acosy, J.=^^smy, (B19) 

we find 

Co= ^dyyJ^{y)J^{y) + ^cos2y, = ^. (B20) 

JO 2 2 ^vr 2tt 



Now we use (|B20|) . (IBTTI) and (IBTSJ) together with equations (IBTsD and (IBT3|) to obtain 



in the limit of small n 

4^i + ilny,. (B21) 

n 2 

Using this together with fl98|) and fllOOp we deduce that 

Bn = e^/". (B22) 



< = -^, (B23) 



Substituting ( 1B22I) in (11091) and taking the limit n — )■ we get 

2mQ 

1 

where we use the fact that Pc/D = 1 as n — )■ 0. Recalling that /i = a/{2Q) we see that ( JB23P 
is equivalent to (1B10I) . 
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